#include "header.h"

int main(int argc, char** argv)
{
  work	  wrk;
  params  ctrl;
  aux	  data;
  int	  Emm, K, Jay, l;
  int 	  nref;

  init_restore(&argc, argv, &ctrl, &wrk);
  init_prep(&ctrl);
  init_full_write(&ctrl, &data);
  restore_x(&ctrl, &wrk);
  //write_ImS(&wrk);
  write_Rew(&wrk);
  write_Imw(&wrk);

  Jay = 256;
  K = lrintq(logq(ctrl.n)/logq(2));
  l = lrintq(logq(Jay)/logq(2));
  Emm = K + 1 - l;
  printf("Refinement parameters:\n");
  printf("K = %d\nJay = %d\nEmm = %d\nl = %d\n", K, Jay, Emm, l);
  
  int skip[Emm];
  __float128  dx[Emm];
  for (int j = 0; j < Emm; j++) {
    skip[j] = (1<<j) ;
    dx[j] = 1.Q*powq(2.,j)/(ctrl.n);
    printf("At level %d of %d: Skip = %d\tdu = %.31Qe\n", j+1, Emm, skip[j], dx[j]);
  }
  nref = 1 + (1<<(l-1))*(Emm+1);
  printf("Total points: %d\n", nref-2);
  __float128 x_pos[nref];
  __float128 XX[nref];
  __float128 YY[nref];
  x_pos[0] = -0.5Q;
  XX[0] = data.XX[0];
  YY[0] = data.YY[0];
  int kk = 0;
  int index = 0;

  for (int mm = 0; mm < Emm; mm++) {
    for (int nn = 0; nn < Jay/4; nn++) {
      kk++;
      index = index + skip[Emm-mm-1];
      x_pos[kk] = x_pos[kk-1] + 1.Q*dx[Emm-mm-1];
      XX[kk] = data.XX[index];
      YY[kk] = data.YY[index];
    }
  }
  for (int nn = Jay/4; nn < Jay/2+1; nn++) {
    kk++; index++;
    x_pos[kk] = x_pos[kk-1] + 1.Q*dx[0];
    XX[kk] = data.XX[index];
    YY[kk] = data.YY[index];
  }
  x_pos[kk] = 0.Q;

  //XX[kk+3] = data.XX[index+3];
  //YY[kk+3] = data.YY[index+3];

  for (int j = 0; j < (nref-1)/2; j++) {
    x_pos[nref-1-j] = -1.Q*x_pos[j];
    XX[nref-1-j] = -1.Q*XX[j];
    YY[nref-1-j] = 1.Q*YY[j];
  }

  FILE    *fh1, *fh2, *fh3, *fh_uni; 
  fh1 = fopen("rawImS.dat","w");
  fh2 = fopen("rawImw.dat","w");
  fh3 = fopen("rawRew.dat","w");
  fh_uni = fopen("rawUni.dat","w");
  fprintf(fh1, "\n");
  fprintf(fh2, "\n");
  fprintf(fh3, "\n");
  for (int j = 1; j < nref-1; j++) {
    fprintf(fh_uni, "%.31Qe\n", 0.5Q*x_pos[j] );
    fprintf(fh1, "%.31Qe\n", tanq(1.Q*pi*x_pos[j]) );
    fprintf(fh2, "%.31Qe\n", XX[j]);
    fprintf(fh3, "%.31Qe\n", YY[j]);
  }
  fclose(fh1);
  fclose(fh2);
  fclose(fh3);
  fclose(fh_uni);
}
